function [celltab,xsim]=plotbeta(alpha,beta,into); 
%function plotbeta(alpha,beta,into); 
% Plots the pdf of an B(alpha,beta) distribution 
%   including mean and variance 
% Legend provides percentiles 
% into.save   =1 save 
% into.fpath  file path 
% into.opath  output path
% Output
% -------
%  celltab   cell, contains percentiles 
%            First two rows contain mean and SD respectively
%  xsim      vector of simulated values used to compute percentiles 
%
% Modified 6/02/03
% ------------------------------------------------------------------
if nargin < 3 | isempty(into) == 1 
    into.save=0; 
end; 

meant=alpha/(alpha+beta); 
vart=alpha*beta/(((alpha+beta)^2)*(alpha+beta+1));
sd=sqrt(vart); 
x=0:0.005:1; 
pdf=betapdf(x,alpha,beta); 
figure; 
plot(x,pdf); 
XLim([0 1.05])
astr=num2str(alpha); 
bstr=num2str(beta); 
mstr=num2str(meant); 
sstr=num2str(sd); 
t1=['Beta(' astr ',' bstr ')   mean : ' mstr ' sd: ' sstr ]; 
title(t1);
xsim=betarnd(alpha,beta,20000,1); 
xsim=sort(xsim); 
med=median(xsim);
per=[0.001 0.01 0.05 0.1 0.5 0.9 0.95 0.99 0.999]; 
perc=num2cell(per');
perc2=betainv(per,alpha,beta); 
perc2=num2cell(perc2'); 
jj=1; tstr=length(perc); 
celltit=cell(tstr,1); 
celltab=cell(tstr,2); 
for jj=1:tstr; 
    celltit{jj}=[perc{jj} '    ' perc2{jj}];
    celltab{jj,1}=perc{jj}; 
    celltab{jj,2}=perc2{jj};
end 
legend('density',celltit{1},celltit{2},celltit{3},celltit{4},celltit{5},celltit{6},celltit{7},celltit{8},celltit{9})
celladd=[{'1'} {mstr}; {'2'} {sstr}]; 
celltab=[celladd;celltab]; 
if into.save==1 
    into.desc=['Beta' astr '_' bstr]; 
    savegraph(into)
end